Important Bird Areas

library(dplyr)
library(here)
library(mapview)
library(tigris)
library(sf)

mapviewOptions(basemaps.color.shuffle = FALSE)

ca.blocks.clean <- readRDS(here("Data", "CABlocks", "CABlocks2024.RDS")) %>%
  filter(WATER_SLIVER == 0 & CA_SLIVER == 0)
cntr_crds <- c(mean(st_coordinates(ca.blocks.clean)[, 1]),
               mean(st_coordinates(ca.blocks.clean)[, 2]))

# Load IBA data
IBAs <- st_read(here("Data", "IBAs", "iba_polygons_public", "iba_polygons_public.shp")) %>%
  st_transform(st_crs(ca.blocks.clean)) %>%
  st_make_valid() %>%
  mutate(IBA_area = st_area(.) %>% units::drop_units())

# Intersect with CA-BBA blocks, count number of IBAs by block, and record thier IDs
iba_intersect <- st_intersection(ca.blocks.clean, IBAs)

iba_intersect <- iba_intersect %>%
  mutate(intersection_area = st_area(geometry)) %>%
  as.data.frame() %>%
  st_drop_geometry() %>%
  select(BLOCKCODE, site_id, intersection_area) %>%
  mutate(intersection_area = units::drop_units(intersection_area)) %>%
  group_by(BLOCKCODE, site_id) %>%
  summarize(intersection_area = sum(intersection_area), .groups = 'drop') %>%
  group_by(BLOCKCODE) %>%
  summarize(
    IBA_COUNT = length(unique(site_id)),
    IBA_ALL_ID = paste(unique(site_id), collapse = ", "),
    IBA_MOST_COVERAGE = site_id[which.max(intersection_area)],
    .groups = 'drop'
  )

# For blocks with multiple IBAs, identify which IBA is the most localized (smallest total area). Also pull names of all IBAs in each block using IBA ids.
iba_intersect$IBA_MOST_LOCALIZED_ID <- NA
iba_intersect$IBA_ALL_NAMES <- NA
for(i in 1:nrow(iba_intersect)){
  rel_ibas <- IBAs %>% st_drop_geometry() %>% filter(site_id %in% as.numeric(unlist(strsplit(iba_intersect$IBA_ALL_ID[i],split=", ",fixed=TRUE))))
  iba_intersect[i, "IBA_MOST_LOCALIZED_ID"] <- rel_ibas %>% slice_min(IBA_area) %>% pull(site_id)
  iba_intersect[i, "IBA_ALL_NAMES"] <- paste(rel_ibas %>% pull(site_name), collapse = "; ")
}

# Join IBA info back to full blocks
ca.blocks.clean <- ca.blocks.clean %>%
  left_join(iba_intersect, by = "BLOCKCODE") %>%
  left_join(IBAs %>% st_drop_geometry() %>% select(site_id, site_name, priority), by = c("IBA_MOST_LOCALIZED_ID" = "site_id")) %>%
  rename("IBA_MOST_LOCALIZED_NAME" = "site_name",
         "IBA_MOST_LOCALIZED_PRIORITY" = "priority")